program main
	!本程序提供了一个实现，该实现调用线性方程组求解函数依次求解算例文件夹中的所有方程组，并输出求解耗时和所有线性方程组求解结果的最大残差
	use ReadData, only: readMatrix, readVector
	USE LAPACK95
	use WriteData, only:outputSolution
	use ComputeResidual, only:computeLargestRelResid
	use CheckData, only:matrixMatchExplanation
	
implicit none
	character(len=100) :: folderName = "StepMatrixRHSAndSolution-3/"!算例的文件夹名
	integer, parameter :: fileNum = 46!算例中的线性方程组个数
	
	integer, parameter ::  nTestRounds = 80 !每个线性方程组求解测试的轮数，此处暂给为1，实际测试过程过程中会取平均求解时间达到稳定时的轮数
	
	!与数据文件名相关的变量和参数，不要修改
	character(len=5) :: matrixSuffix = ".matr"
	character(len=4) :: rhsSuffix = ".rhs"
	character(len=100) :: formatDescription
	character(len=100) :: fileNamePrefixM, fileNamePrefixN, fileNamePrefix
	integer, parameter :: fileNamePrefixLength = 5
	
	real(kind=8), dimension(:,:), pointer :: matrix !线性方程组的矩阵
	real(kind=8), dimension(:), pointer :: rhs !线性方程组的右端项
	
	integer :: i,j,k
	
	!用于计时的相关参数
	integer :: count1, count2, countRate 
	integer :: totalCount = 0
	
	!与localMatrix和localRhs相关的参数，定义这些变量是因为gemv会修改传入的矩阵和右端项
	integer :: matrixSize
	real(kind=8), dimension(:,:), allocatable :: localMatrix
	real(kind=8), dimension(:), allocatable :: localRhs
	
	write(formatDescription,'(I1)') fileNamePrefixLength
	formatDescription = '(I'//TRIM(ADJUSTL(formatDescription))//')'
	
	!检查该算例数据中的数据是否满足赛题描述中的分块矩阵前导子矩阵保持不变的性质
	!此部分主要供选手参考矩阵性质，在实际测试选手的成果时会注释掉
	!***************************************************************************************
	
	do i = 1,fileNum - 1
		write(fileNamePrefixM,formatDescription) i
		
		do k = 1, fileNamePrefixLength - LEN(TRIM(ADJUSTL(fileNamePrefixM)))
			fileNamePrefixM  = '0'//ADJUSTL(TRIM(fileNamePrefixM))
		end do
 
		do j = i+1, fileNum
			
			write(fileNamePrefixN,formatDescription) j
			do k = 1, fileNamePrefixLength - LEN(TRIM(ADJUSTL(fileNamePrefixN)))
				fileNamePrefixN  = '0'//ADJUSTL(TRIM(fileNamePrefixN))
			end do
			
			if (.not. matrixMatchExplanation(TRIM(ADJUSTL(folderName)) // TRIM(ADJUSTL(fileNamePrefixM)), TRIM(ADJUSTL(folderName)) // TRIM(ADJUSTL(fileNamePrefixN)), matrixSuffix) ) then
				print *, "Data not matching explanation"
				pause
				stop
			end if
		end do
	end do
	print *,"Data matching explanation"
	pause
	
	!***************************************************************************************
	
	
	do i = 1, fileNum
		write(fileNamePrefix,formatDescription) i
		do j = 1, fileNamePrefixLength - LEN(TRIM(ADJUSTL(fileNamePrefix)))
			fileNamePrefix  = '0'//ADJUSTL(TRIM(fileNamePrefix))
		end do
		
		call readMatrix(TRIM(ADJUSTL(folderName)) // TRIM(ADJUSTL(fileNamePrefix)),matrixSuffix,matrix)
		call readVector(TRIM(ADJUSTL(folderName)) // TRIM(ADJUSTL(fileNamePrefix)),rhsSuffix,rhs)
		
		matrixSize = size(rhs,1)
		allocate(localMatrix(matrixSize,matrixSize), localRhs(matrixSize))
		
		do j = 1,nTestRounds
			localMatrix(:,:) = matrix(:,:)
			localRhs(:) = rhs(:)
			call SYSTEM_CLOCK(count1)!本实现采用fortran内置的SYSTEM_CLOCK进行计时，会使用统一的计时子程序对所有选手的作品进行评估
			call gesv(localMatrix,localRhs) !本实现调用了mkl 中的gemv来对线性方程组进行求解，后续此部分将会替换成参赛者提交的程序。注意：由于参赛者可能会用到矩阵逐步扩大但前导子矩阵不变的性质，因此程序接口不用完全与gemv一样。
			call SYSTEM_CLOCK(count2)
			totalCount = totalCount + count2 - count1
		end do
		
		call outputSolution(localRhs,folderName,fileNamePrefix)
		
		deallocate (matrix,rhs,localMatrix,localRhs)
		
	end do
	
	call SYSTEM_CLOCK(count_rate=countRate)
	print *, 'Time taken: ', real(totalCount) / real(countRate), ' seconds'
	
	print *, 'Largest relative residual: ', computeLargestRelResid(folderName,fileNum)
	
	pause


end program main
